home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cbsp_aux.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  17.7 KB  |  568 lines
  1. /******************************************************************************
  2. * CBsp-Aux.c - Bspline curve auxilary routines.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /******************************************************************************
  13. * Given a bspline curve - subdivide it into two at the given parametric value.*
  14. * Returns pointer to first curve in a list of two curves (subdivided ones).   *
  15. * The subdivision is achieved by inserting (order-1) knot at the given param. *
  16. * value t and spliting the control polygon and knot vector at that point.     *
  17. ******************************************************************************/
  18. CagdCrvStruct *BspCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
  19. {
  20.     CagdBType
  21.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  22.     int i, j,
  23.     k = Crv -> Order,
  24.     Len = Crv -> Length,
  25.         KVLen = k + Len,
  26.     Index1 = BspKnotLastIndexL(Crv -> KnotVector, KVLen, t),
  27.     Index2 = BspKnotFirstIndexG(Crv -> KnotVector, KVLen, t),
  28.     Mult = k - 1 - (Index2 - Index1 - 1),
  29.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  30.     CagdRType TMin, TMax, *LOnePts, *ROnePts, *OnePts, *NewKV,
  31.     **Pts, **LPts, **RPts;
  32.     CagdCrvStruct *LCrv, *RCrv;
  33.     BspKnotAlphaCoeffType *A;
  34.  
  35.     CagdCrvDomain(Crv, &TMin, &TMax);
  36.     if (t < TMin || APX_EQ(t, TMin) ||
  37.     t > TMax || APX_EQ(t, TMax))
  38.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  39.  
  40.     LCrv = BspCrvNew(Index1 + 1, k, Crv -> PType);
  41.     RCrv = BspCrvNew(Len - Index2 + k, k, Crv -> PType);
  42.  
  43.     /* Update the new knot vectors. */
  44.     CAGD_GEN_COPY(LCrv -> KnotVector,
  45.          Crv -> KnotVector,
  46.          sizeof(CagdRType) * (Index1 + 1));
  47.     /* Close the knot vector with multiplicity Order: */
  48.     for (j = Index1 + 1; j <= Index1 + k; j++)
  49.     LCrv -> KnotVector[j] = t;
  50.     CAGD_GEN_COPY(&RCrv -> KnotVector[k],
  51.           &Crv -> KnotVector[Index2],
  52.           sizeof(CagdRType) * (Len + k - Index2));
  53.     /* Make sure knot vector starts with multiplicity Order: */
  54.     for (j = 0; j < k; j++)
  55.     RCrv -> KnotVector[j] = t;
  56.  
  57.     /* Now handle the control polygon refinement. */
  58.     Pts = Crv -> Points,
  59.     LPts = LCrv -> Points,
  60.     RPts = RCrv -> Points;
  61.  
  62.     if (Mult > 0) {
  63.     NewKV = IritMalloc(sizeof(CagdRType) * Mult);
  64.     for (i = 0; i < Mult; i++)
  65.         NewKV[i] = t;
  66.     A = BspKnotEvalAlphaCoefMerge(k, Crv -> KnotVector, Len, NewKV,
  67.                       Mult);
  68.     IritFree((VoidPtr) NewKV);
  69.     }
  70.     else {
  71.     A = BspKnotEvalAlphaCoef(k, Crv -> KnotVector, Len,
  72.                     Crv -> KnotVector, Len);
  73.     }
  74.  
  75.     /* Note that Mult can be negative in cases where original       */
  76.     /* multiplicity was order or more and we need to compensate     */
  77.     /* here, since Alpha matrix will be just a unit matrix then.    */
  78.     Mult = Mult >= 0 ? 0 : -Mult;
  79.  
  80.     /* Blend Crv into LCrv. */
  81.     for (j = IsNotRational; j <= MaxCoord; j++) {
  82.     LOnePts = LPts[j];
  83.     OnePts = Pts[j];
  84.     for (i = 0; i < LCrv -> Length; i++, LOnePts++)
  85.         CAGD_ALPHA_BLEND( A, i, OnePts, LOnePts );
  86.     }
  87.  
  88.     /* Blend Crv into RCrv. */
  89.     for (j = IsNotRational; j <= MaxCoord; j++) {
  90.     ROnePts = RPts[j];
  91.     OnePts = Pts[j];
  92.     for (i = LCrv -> Length - 1 + Mult;
  93.          i < LCrv -> Length + RCrv -> Length - 1 + Mult;
  94.          i++, ROnePts++)
  95.         CAGD_ALPHA_BLEND( A, i, OnePts, ROnePts );
  96.     }
  97.  
  98.     BspKnotFreeAlphaCoef(A);
  99.  
  100.     BspKnotMakeRobustKV(RCrv -> KnotVector,
  101.             RCrv -> Order + RCrv -> Length);
  102.     BspKnotMakeRobustKV(LCrv -> KnotVector,
  103.             LCrv -> Order + LCrv -> Length);
  104.  
  105.     LCrv -> Pnext = RCrv;
  106.  
  107.     return LCrv;
  108. }
  109.  
  110. /******************************************************************************
  111. *  Insert n knot all with the value t. In no case will the multiplicity of    *
  112. * knot be greater or equal to the curve order.                      *
  113. ******************************************************************************/
  114. CagdCrvStruct *BspCrvKnotInsertNSame(CagdCrvStruct *Crv, CagdRType t, int n)
  115. {
  116.     int i,
  117.     CrntMult = BspKnotFindMult(Crv -> KnotVector, Crv -> Order,
  118.                                Crv -> Length, t),
  119.     Mult = MIN(n, Crv -> Order - CrntMult - 1);
  120.     CagdCrvStruct *RefinedCrv;
  121.  
  122.     if (Mult > 0) {
  123.     CagdRType
  124.         *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);
  125.  
  126.     for (i = 0; i < Mult; i++)
  127.         NewKV[i] = t;
  128.  
  129.     RefinedCrv = BspCrvKnotInsertNDiff(Crv, FALSE, NewKV, Mult);
  130.  
  131.     IritFree((VoidPtr) NewKV);
  132.     }
  133.     else {
  134.     RefinedCrv = CagdCrvCopy(Crv);
  135.     }
  136.  
  137.     return RefinedCrv;
  138. }
  139.  
  140. /******************************************************************************
  141. *  Insert n knot with different values as defined by t. If however Replace is *
  142. * TRUE, the knot are simply replacing the current ones.                  *
  143. ******************************************************************************/
  144. CagdCrvStruct *BspCrvKnotInsertNDiff(CagdCrvStruct *Crv, CagdBType Replace,
  145.                                CagdRType *t, int n)
  146. {
  147.     CagdBType
  148.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  149.     CagdRType
  150.     *KnotVector = Crv -> KnotVector;
  151.     int Order = Crv -> Order,
  152.     Length = Crv -> Length,
  153.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  154.     CagdCrvStruct *RefinedCrv;
  155.  
  156.     if (Replace) {
  157.     int i;
  158.  
  159.     if (Order + Length != n)
  160.         FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
  161.     for (i = 1; i < n; i++)
  162.         if (t[i] < t[i - 1])
  163.         FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
  164.  
  165.         RefinedCrv = CagdCrvCopy(Crv);
  166.     for (i = 0; i < n; i++)
  167.         RefinedCrv -> KnotVector[i] = *t++;
  168.     }
  169.     else if (n == 0) {
  170.     RefinedCrv = CagdCrvCopy(Crv);
  171.     }
  172.     else {
  173.     int i, j, LengthKVt;
  174.     BspKnotAlphaCoeffType *A;
  175.     CagdRType *MergedKVt;
  176.  
  177.     /* Compute the Alpha refinement matrix. */
  178.     MergedKVt = BspKnotMergeTwo(KnotVector, Length + Order,
  179.                     t, n, 0, &LengthKVt);
  180.     A = BspKnotEvalAlphaCoef(Order, KnotVector, Length,
  181.                  MergedKVt, LengthKVt - Order);
  182.  
  183.     RefinedCrv = BspCrvNew(Crv -> Length + n, Order, Crv -> PType);
  184.  
  185.     /* Update the knot vector. */
  186.     IritFree((VoidPtr) RefinedCrv -> KnotVector);
  187.     RefinedCrv -> KnotVector = MergedKVt;
  188.  
  189.     /* Update the control mesh */
  190.     for (j = IsNotRational; j <= MaxCoord; j++) {
  191.         CagdRType
  192.         *ROnePts = RefinedCrv -> Points[j],
  193.         *OnePts = Crv -> Points[j];
  194.         for (i = 0; i < RefinedCrv -> Length; i++, ROnePts++)
  195.         CAGD_ALPHA_BLEND( A, i, OnePts, ROnePts );
  196.     }
  197.     BspKnotFreeAlphaCoef(A);
  198.     }
  199.  
  200.     BspKnotMakeRobustKV(RefinedCrv -> KnotVector,
  201.             RefinedCrv -> Order + RefinedCrv -> Length);
  202.  
  203.     return RefinedCrv;
  204. }
  205.  
  206. /******************************************************************************
  207. * Returns a new curve, identical to the original but with order N.          *
  208. *   Degree raise by multiplying by a constant 1 curve of order                  *
  209. * (NewOrder - Order).                                  *
  210. ******************************************************************************/
  211. CagdCrvStruct *BspCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
  212. {
  213.     int i, j, RaisedOrder,
  214.     Order = Crv -> Order,
  215.     Length = Crv -> Length,
  216.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType),
  217.     KvLen1 = Order + Length - 1;
  218.     CagdCrvStruct *RaisedCrv, *UnitCrv;
  219.     CagdRType
  220.     *Kv = Crv -> KnotVector;
  221.  
  222.     if (NewOrder < Order) {
  223.     FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
  224.     return NULL;
  225.     }
  226.     RaisedOrder = NewOrder - Order + 1;
  227.  
  228.     UnitCrv = BspCrvNew(RaisedOrder, RaisedOrder,
  229.             CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  230.     for (i = 0; i < RaisedOrder * 2; i++)
  231.     UnitCrv -> KnotVector[i] = i >= RaisedOrder ? Kv[KvLen1] : Kv[0];
  232.     for (i = 1; i <= MaxCoord; i++)
  233.     for (j = 0; j < RaisedOrder; j++)
  234.         UnitCrv -> Points[i][j] = 1.0;
  235.  
  236.     RaisedCrv = BspCrvMult(Crv, UnitCrv);
  237.  
  238.     CagdCrvFree(UnitCrv);
  239.  
  240.     return RaisedCrv;
  241. }
  242.  
  243. /******************************************************************************
  244. * Return a new curve, identical to the original but with one degree higher    *
  245. ******************************************************************************/
  246. CagdCrvStruct *BspCrvDegreeRaise(CagdCrvStruct *Crv)
  247. {
  248.     CagdBType
  249.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  250.     int i, i2, j, RaisedLen,
  251.     Order = Crv -> Order,
  252.     Length = Crv -> Length,
  253.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  254.     CagdCrvStruct *RaisedCrv;
  255.  
  256.     if (Order > 2)
  257.     return BspCrvDegreeRaiseN(Crv, Order + 1);
  258.  
  259.     /* If curve is linear, degree raising means basically to increase the    */
  260.     /* knot multiplicity of each segment by one and add a middle point for   */
  261.     /* each such segment.                             */
  262.     RaisedLen = Length * 2 - 1;
  263.     RaisedCrv = BspCrvNew(RaisedLen, Order + 1, Crv -> PType);
  264.  
  265.     /* Update the control polygon; */
  266.     for (j = IsNotRational; j <= MaxCoord; j++)             /* First point. */
  267.     RaisedCrv -> Points[j][0] = Crv -> Points[j][0];
  268.  
  269.     for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
  270.     for (j = IsNotRational; j <= MaxCoord; j++) {
  271.         RaisedCrv -> Points[j][i2] =
  272.         Crv -> Points[j][i-1] * 0.5 + Crv -> Points[j][i] * 0.5;
  273.         RaisedCrv -> Points[j][i2 + 1] = Crv -> Points[j][i];
  274.     }
  275.  
  276.     /* Update the knot vector. */
  277.     for (i = 0; i < 3; i++)
  278.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[0];
  279.  
  280.     for (i = 2, j = 3; i < Length; i++, j += 2)
  281.     RaisedCrv -> KnotVector[j] = RaisedCrv -> KnotVector[j + 1] = 
  282.         Crv -> KnotVector[i];
  283.  
  284.     for (i = j; i < j + 3; i++)
  285.     RaisedCrv -> KnotVector[i] = Crv -> KnotVector[Length];
  286.  
  287.     return RaisedCrv;
  288. }
  289.  
  290. /******************************************************************************
  291. * Return a normalized vector, equal to the tangent to Crv at parameter t.     *
  292. * Algorithm: insert (order - 1) knots and return control polygon tangent.     *
  293. ******************************************************************************/
  294. CagdVecStruct *BspCrvTangent(CagdCrvStruct *Crv, CagdRType t)
  295. {
  296.     static CagdVecStruct P2;
  297.     CagdVecStruct P1, *T;
  298.     CagdRType TMin, TMax;
  299.     int k = Crv -> Order,
  300.     Len = Crv -> Length,
  301.     Index = BspKnotLastIndexL(Crv -> KnotVector, k + Len, t);
  302.     CagdPointType
  303.     PType = Crv -> PType;
  304.     CagdCrvStruct *RefinedCrv;
  305.  
  306.     if (!BspKnotParamInDomain(Crv -> KnotVector, Len, k, t))
  307.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  308.  
  309.     BspCrvDomain(Crv, &TMin, &TMax);
  310.  
  311.     if (APX_EQ(t, TMin)) {
  312.     /* Use Crv starting tangent direction. */
  313.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, PType);
  314.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, PType);
  315.     }
  316.     else if (APX_EQ(t, TMax)) {
  317.     /* Use Crv ending tangent direction. */
  318.     CagdCoerceToE3(P1.Vec, Crv -> Points, Len - 2, PType);
  319.     CagdCoerceToE3(P2.Vec, Crv -> Points, Len - 1, PType);
  320.     }
  321.     else {
  322.     RefinedCrv = BspCrvKnotInsertNSame(Crv, t, k - 1);
  323.  
  324.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  325.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  326.  
  327.     CagdCrvFree(RefinedCrv);
  328.     }
  329.  
  330.     CAGD_SUB_VECTOR(P2, P1);
  331.  
  332.     if (CAGD_LEN_VECTOR(P2) < SQR(EPSILON)) {
  333.     if (AttrGetIntAttrib(Crv -> Attr, "_tan") != TRUE) {
  334.         /* Try to move a little. This location has zero speed. However,  */
  335.         /* do it only once since we can be here forever. The "_tan"      */
  336.         /* attribute guarantee we will try to move EPSILON only once.    */
  337.         AttrSetIntAttrib(&Crv -> Attr, "_tan", TRUE);
  338.  
  339.         T = BspCrvTangent(Crv, t - TMin < TMax - t ? t + EPSILON
  340.                                : t - EPSILON);
  341.  
  342.         AttrFreeOneAttribute(&Crv -> Attr, "_tan");
  343.  
  344.         return T;
  345.     }
  346.     else {
  347.         /* A zero length vector signals failure to compute tangent. */
  348.         return &P2;
  349.     }
  350.     }
  351.     else {
  352.     CAGD_NORMALIZE_VECTOR(P2);            /* Normalize the vector. */
  353.  
  354.     return &P2;
  355.     }
  356. }
  357.  
  358. /******************************************************************************
  359. * Return a normalized vector, equal to the binormal to Crv at parameter t.    *
  360. * Algorithm: insert (order - 1) knots and using 3 consecutive control points  *
  361. * at the refined location (p1, p2, p3), compute to binormal to be the cross   *
  362. * product of the two vectors (p1 - p2) and (p2 - p3).                  *
  363. ******************************************************************************/
  364. CagdVecStruct *BspCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
  365. {
  366.     static CagdVecStruct P3;
  367.     CagdVecStruct P1, P2;
  368.     CagdRType TMin, TMax;
  369.     int k = Crv -> Order,
  370.     Len = Crv -> Length,
  371.     Index = BspKnotLastIndexL(Crv -> KnotVector, k + Len, t);
  372.     CagdPointType
  373.     PType = Crv -> PType;
  374.     CagdCrvStruct *RefinedCrv;
  375.  
  376.     if (!BspKnotParamInDomain(Crv -> KnotVector, Len, k, t))
  377.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  378.  
  379.     /* Can not compute for linear curves. */
  380.     if (k <= 2)
  381.     return NULL;
  382.  
  383.     /* For planar curves, B is trivially the Z axis. */
  384.     if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
  385.     P3.Vec[0] = P3.Vec[1] = 0.0;
  386.     P3.Vec[2] = 1.0;
  387.     return &P3;
  388.     }
  389.  
  390.     BspCrvDomain(Crv, &TMin, &TMax);
  391.  
  392.     if (APX_EQ(t, TMin)) {
  393.     /* Use Crv starting tangent direction. */
  394.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, PType);
  395.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, PType);
  396.     CagdCoerceToE3(P3.Vec, Crv -> Points, 2, PType);
  397.     }
  398.     else if (APX_EQ(t, TMax)) {
  399.     /* Use Crv ending tangent direction. */
  400.     CagdCoerceToE3(P1.Vec, Crv -> Points, Len - 3, PType);
  401.     CagdCoerceToE3(P2.Vec, Crv -> Points, Len - 2, PType);
  402.     CagdCoerceToE3(P3.Vec, Crv -> Points, Len - 1, PType);
  403.     }
  404.     else {
  405.     RefinedCrv = BspCrvKnotInsertNSame(Crv, t, k - 1);
  406.  
  407.     CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
  408.     CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
  409.     CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, Index + 2, PType);
  410.  
  411.     CagdCrvFree(RefinedCrv);
  412.     }
  413.  
  414.     CAGD_SUB_VECTOR(P1, P2);
  415.     CAGD_SUB_VECTOR(P2, P3);
  416.  
  417.     CROSS_PROD(P3.Vec, P1.Vec, P2.Vec);
  418.  
  419.     if ((t = CAGD_LEN_VECTOR(P3)) < EPSILON)
  420.     return NULL;
  421.     else
  422.     CAGD_DIV_VECTOR(P3, t);                /* Normalize the vector. */
  423.  
  424.     return &P3;
  425. }
  426.  
  427. /******************************************************************************
  428. * Return a normalized vector, equal to the normal to Crv at parameter t.      *
  429. * Algorithm: returns the cross product of the curve tangent and binormal.     *
  430. ******************************************************************************/
  431. CagdVecStruct *BspCrvNormal(CagdCrvStruct *Crv, CagdRType t)
  432. {
  433.     static CagdVecStruct N, *T, *B;
  434.  
  435.     T = BspCrvTangent(Crv, t);
  436.     B = BspCrvBiNormal(Crv, t);
  437.  
  438.     if (T == NULL || B == NULL)
  439.     return NULL;
  440.  
  441.     CROSS_PROD(N.Vec, T -> Vec, B -> Vec);
  442.  
  443.     CAGD_NORMALIZE_VECTOR(N);                /* Normalize the vector. */
  444.  
  445.     return &N;
  446. }
  447.  
  448. /******************************************************************************
  449. * Return a new curve, equal to the derived curve.                  *
  450. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
  451. * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.   *
  452. ******************************************************************************/
  453. CagdCrvStruct *BspCrvDerive(CagdCrvStruct *Crv)
  454. {
  455.     CagdBType
  456.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  457.     int i, j,
  458.     k = Crv -> Order,
  459.     Len = Crv -> Length,
  460.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  461.     CagdRType
  462.     *Kv = Crv -> KnotVector;
  463.     CagdCrvStruct *DerivedCrv;
  464.  
  465.     if (!IsNotRational)
  466.     return BspCrvDeriveRational(Crv);
  467.     else if (k < 2)
  468.     FATAL_ERROR(CAGD_ERR_LIN_NO_SUPPORT);
  469.  
  470.     DerivedCrv = BspCrvNew(Len - 1, k - 1, Crv -> PType);
  471.  
  472.     for (i = 0; i < Len - 1; i++) {
  473.     CagdRType
  474.         Denom = Kv[i + k] - Kv[i + 1];
  475.  
  476.     for (j = IsNotRational; j <= MaxCoord; j++)
  477.         DerivedCrv -> Points[j][i] = (k - 1) *
  478.         (Crv -> Points[j][i + 1] - Crv -> Points[j][i]) / Denom;
  479.     }
  480.  
  481.     CAGD_GEN_COPY(DerivedCrv -> KnotVector,
  482.           &Crv -> KnotVector[1],
  483.           sizeof(CagdRType) * (k + Len - 2));
  484.  
  485.     return DerivedCrv;
  486. }
  487.  
  488. /******************************************************************************
  489. * Return a new Bspline curve, equal to the integral of the given Bspline crv. *
  490. * The given Bspline curve should be nonrational.                  *
  491. *                                          *
  492. *           l           l           l       l+1                  *
  493. *   /         /-           -      /       -   P    -                  *
  494. *  |        | \        n       \     |  n       \    i   \              n+1          *
  495. *  | C(t) = | / P  B (t) = / P   | B (t) = / -----  / ( t   - t ) B   (t) =   *
  496. * /        /  -     i  i       -  i /   i       - n + 1  -    j+n   j   j          *
  497. *            i=0          i=0             i=0     j=i+1                  *
  498. *                                          *
  499. *        l+1 j-1                                  *
  500. *         -   -                                      *
  501. *     1   \   \            n+1                          *
  502. * = ----- /   / P  ( t   - t ) B   (t)                          *
  503. *   n + 1 -   -  i    j+n   j    j                          *
  504. *        j=1 i=0                                  *
  505. *                                          *
  506. ******************************************************************************/
  507. CagdCrvStruct *BspCrvIntegrate(CagdCrvStruct *Crv)
  508. {
  509.     int i, j, k,
  510.     n = Crv -> Order,
  511.     Len = Crv -> Length,
  512.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  513.     CagdCrvStruct *IntCrv;
  514.     CagdRType
  515.     *Kv = Crv -> KnotVector;
  516.  
  517.     if (CAGD_IS_RATIONAL_CRV(Crv))
  518.     FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT);
  519.  
  520.     IntCrv = BspCrvNew(Len + 1, n + 1, Crv -> PType);
  521.  
  522.     /* Copy the knot vector and duplicate the two end knots. */
  523.     CAGD_GEN_COPY(&IntCrv -> KnotVector[1], Kv, sizeof(CagdRType) * (Len + n));
  524.     IntCrv -> KnotVector[0] = Kv[0];
  525.     IntCrv -> KnotVector[Len + n + 1] = Kv[Len + n - 1];
  526.     Kv = IntCrv -> KnotVector;
  527.  
  528.     for (k = 1; k <= MaxCoord; k++) {
  529.     CagdRType
  530.         *Points = Crv -> Points[k],
  531.         *IntPoints = IntCrv -> Points[k];
  532.  
  533.     for (j = 0; j < Len + 1; j++) {
  534.         IntPoints[j] = 0.0;
  535.         for (i = 0; i < j; i++)
  536.             IntPoints[j] += Points[i] * (Kv[i + n + 1] - Kv[i + 1]);
  537.         IntPoints[j] /= n;
  538.     }
  539.     }
  540.  
  541.     return IntCrv;
  542. }
  543.  
  544. /******************************************************************************
  545. * Return a new linear Bspline curve constructed from the given polyline.      *
  546. ******************************************************************************/
  547. CagdCrvStruct *CnvrtPolyline2LinBsplineCrv(CagdPolylineStruct *Poly)
  548. {
  549.     int i,
  550.     Length = Poly -> Length;
  551.     CagdCrvStruct
  552.     *Crv = BspCrvNew(Length, 2, CAGD_PT_E3_TYPE);
  553.     CagdRType
  554.     **Points = Crv -> Points;
  555.     CagdPtStruct
  556.     *Pts = Poly -> Polyline;
  557.  
  558.     BspKnotUniformOpen(Length, 2, Crv -> KnotVector);
  559.  
  560.     for (i = 0; i < Length; i++, Pts++) {
  561.     Points[1][i] = Pts -> Pt[0];
  562.     Points[2][i] = Pts -> Pt[1];
  563.     Points[3][i] = Pts -> Pt[2];
  564.     }
  565.  
  566.     return Crv;
  567. }
  568.